/*-*-mode: 

!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/* fun-pw86x.c:

   Implemented by Dave Wilson (david.wilson@latrobe.edu.au), Jun 2005

   Automatically generated code implementing PW86X functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional pw86xFunctional;" to 'functionals.h'
    2. add "&pw86xFunctional," to 'functionals.c'
    3. add "fun-pw86x.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
rho:  rhoa + rhob;
grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
zeta: (rhoa-rhob)/(rhoa+rhob);

kfa: (6*%PI*%PI*rhoa)^(1/3);
kfb: (6*%PI*%PI*rhob)^(1/3);
exa: -3*kfa/(4*%PI);
exb: -3*kfb/(4*%PI);

sa:  grada/(2*kfa*rhoa);
sb:  gradb/(2*kfb*rhob);
Fa:  (1+1.296*sa^2+14.0*sa^4+0.2*sa^6)^(1/15);
Fb:  (1+1.296*sb^2+14.0*sb^4+0.2*sb^6)^(1/15);

Exa: 2*rhoa*exa*Fa;
Exb: 2*rhob*exb*Fb;

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb)/2;

    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"
 
/* INTERFACE PART */
static integer pw86x_isgga(void) { return 1; } /* FIXME: detect! */
static integer pw86x_read(const char *conf_line);
static real pw86x_energy(const FunDensProp* dp);
static void pw86x_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void pw86x_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void pw86x_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void pw86x_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional PW86xFunctional = {
  "PW86x",       /* name */
  pw86x_isgga,   /* gga-corrected */
   1,
  pw86x_read,
  NULL,
  pw86x_energy,
  pw86x_first,
  pw86x_second,
  pw86x_third,
  pw86x_fourth
};
 
/* IMPLEMENTATION PART */
static integer
pw86x_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
pw86x_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4, t5, t6, t7;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/pow(3.141592653589793,4.0);
    t4 = 1/t1;
    t5 = 1/pow(3.141592653589793,2.666666666666667);
    t6 = 1/pow(6.0,0.666666666666667);
    t7 = 1/pow(3.141592653589793,1.333333333333333);

   /* code */
    res = 0.5*(-1.5*t1*t2*pow(rhob,1.333333333333333)*pow(8.680555555555556E-5*
        t3*pow(gradb,6.0)/pow(rhob,8.0)+0.145833333333333*t4*t5*pow(gradb,
        4.0)/pow(rhob,5.333333333333333)+0.324*t6*t7*pow(gradb,2.0)/
        pow(rhob,2.666666666666667)+1.0,0.066666666666667)-1.5*t1*
        t2*pow(rhoa,1.333333333333333)*pow(8.680555555555556E-5*t3*
        pow(grada,6.0)/pow(rhoa,8.0)+0.145833333333333*t4*t5*pow(grada,
        4.0)/pow(rhoa,5.333333333333333)+0.324*t6*t7*pow(grada,2.0)/
        pow(rhoa,2.666666666666667)+1.0,0.066666666666667));

    return res;
}

static void
pw86x_first_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/pow(3.141592653589793,4.0);
    t4 = pow(grada,6.0);
    t5 = 1/pow(rhoa,8.0);
    t6 = 1/t1;
    t7 = 1/pow(3.141592653589793,2.666666666666667);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = 1/pow(6.0,0.666666666666667);
    t11 = 1/pow(3.141592653589793,1.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,2.666666666666667);
    t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
        t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
    t15 = 1/pow(t14,0.933333333333333);
    t16 = pow(rhoa,1.333333333333333);

   /* code */
    res[0] = 0.5*(-0.1*t1*t15*t16*t2*(-6.944444444444445E-4*
        t3*t4/pow(rhoa,9.0)-0.777777777777778*t6*t7*t8/pow(rhoa,6.333333333333333)-
        0.864*t10*t11*t12/pow(rhoa,3.666666666666667))-2.0*t1*pow(t14,
        0.066666666666667)*t2*pow(rhoa,0.333333333333333));
    res[1] = -0.05*t1*t15*t16*t2*(5.208333333333333E-4*t3*
        t5*pow(grada,5.0)+0.583333333333333*t6*t7*t9*pow(grada,3.0)+
        0.648*t10*t11*grada*t13);
}

static void
pw86x_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    pw86x_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw86x_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
pw86x_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/pow(3.141592653589793,4.0);
    t4 = pow(grada,6.0);
    t5 = 1/pow(rhoa,8.0);
    t6 = 1/t1;
    t7 = 1/pow(3.141592653589793,2.666666666666667);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = 1/pow(6.0,0.666666666666667);
    t11 = 1/pow(3.141592653589793,1.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,2.666666666666667);
    t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
        t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
    t15 = pow(t14,0.066666666666667);
    t16 = pow(rhoa,0.333333333333333);
    t17 = 1/pow(rhoa,9.0);
    t18 = 1/pow(rhoa,6.333333333333333);
    t19 = 1/pow(rhoa,3.666666666666667);
    t20 = -0.864*t10*t11*t12*t19-0.777777777777778*t6*t7*
        t8*t18-6.944444444444445E-4*t3*t4*t17;
    t21 = 1/pow(t14,0.933333333333333);
    t22 = pow(rhoa,1.333333333333333);
    t23 = pow(grada,5.0);
    t24 = pow(grada,3.0);
    t25 = 0.648*t10*t11*grada*t13+0.583333333333333*t6*t7*
        t24*t9+5.208333333333333E-4*t3*t23*t5;
    t26 = 1/pow(t14,1.933333333333333);

   /* code */
    res[0] = 0.5*(-0.1*t1*t2*t20*t21*t22-2.0*t1*t15*t16*t2);
    res[1] = -0.05*t1*t2*t21*t22*t25;
    res[2] = 0.5*(-0.1*t1*t2*t21*t22*(0.00625*t3*t4/pow(rhoa,
        10.0)+4.925925925925925*t6*t7*t8/pow(rhoa,7.333333333333333)+
        3.168*t10*t11*t12/pow(rhoa,4.666666666666667))-0.666666666666667*
        t1*t15*t2/pow(rhoa,0.666666666666667)+0.093333333333333*t1*
        t2*pow(t20,2.0)*t22*t26-0.266666666666667*t1*t16*t2*t20*t21);
    res[3] = 0.5*(0.093333333333333*t1*t2*t20*t22*t25*t26-
        0.133333333333333*t1*t16*t2*t21*t25-0.1*t1*(-1.728*t10*t11*
        grada*t19-3.111111111111111*t6*t7*t24*t18-0.004166666666667*
        t3*t23*t17)*t2*t21*t22);
    res[4] = 0.046666666666667*t1*t2*t22*pow(t25,2.0)*t26-
        0.05*t1*(0.648*t10*t11*t13+1.75*t6*t7*t12*t9+0.002604166666667*
        t3*t8*t5)*t2*t21*t22;

}

static void
pw86x_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    pw86x_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw86x_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static  void
pw86x_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/pow(3.141592653589793,4.0);
    t4 = pow(grada,6.0);
    t5 = 1/pow(rhoa,8.0);
    t6 = 1/t1;
    t7 = 1/pow(3.141592653589793,2.666666666666667);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = 1/pow(6.0,0.666666666666667);
    t11 = 1/pow(3.141592653589793,1.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,2.666666666666667);
    t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
        t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
    t15 = pow(t14,0.066666666666667);
    t16 = pow(rhoa,0.333333333333333);
    t17 = 1/pow(rhoa,9.0);
    t18 = 1/pow(rhoa,6.333333333333333);
    t19 = 1/pow(rhoa,3.666666666666667);
    t20 = -0.864*t10*t11*t12*t19-0.777777777777778*t6*t7*
        t8*t18-6.944444444444445E-4*t3*t4*t17;
    t21 = 1/pow(t14,0.933333333333333);
    t22 = pow(rhoa,1.333333333333333);
    t23 = pow(grada,5.0);
    t24 = pow(grada,3.0);
    t25 = 0.648*t10*t11*grada*t13+0.583333333333333*t6*t7*
        t24*t9+5.208333333333333E-4*t3*t23*t5;
    t26 = 1/pow(rhoa,0.666666666666667);
    t27 = pow(t20,2.0);
    t28 = 1/pow(t14,1.933333333333333);
    t29 = 1/pow(rhoa,10.0);
    t30 = 1/pow(rhoa,7.333333333333333);
    t31 = 1/pow(rhoa,4.666666666666667);
    t32 = 3.168*t10*t11*t12*t31+4.925925925925925*t6*t7*t8*
        t30+0.00625*t3*t4*t29;
    t33 = -1.728*t10*t11*grada*t19-3.111111111111111*t6*t7*
        t24*t18-0.004166666666667*t3*t23*t17;
    t34 = pow(t25,2.0);
    t35 = 0.648*t10*t11*t13+1.75*t6*t7*t12*t9+0.002604166666667*
        t3*t8*t5;
    t36 = 1/pow(t14,2.933333333333333);

   /* code */
    res[0] = 0.5*(-0.1*t1*t2*t20*t21*t22-2.0*t1*t15*t16*t2);
    res[1] = -0.05*t1*t2*t21*t22*t25;
    res[2] = 0.5*(-0.1*t1*t2*t21*t22*t32+0.093333333333333*
        t1*t2*t22*t27*t28-0.666666666666667*t1*t15*t2*t26-0.266666666666667*
        t1*t16*t2*t20*t21);
    res[3] = 0.5*(-0.1*t1*t2*t21*t22*t33+0.093333333333333*
        t1*t2*t20*t22*t25*t28-0.133333333333333*t1*t16*t2*t21*t25);
    res[4] = 0.046666666666667*t1*t2*t22*t28*t34-0.05*t1*
        t2*t21*t22*t35;
    res[5] = 0.5*(-0.1*t1*t2*t21*t22*(-0.0625*t3*t4/pow(rhoa,
        11.0)-36.12345679012345*t6*t7*t8/pow(rhoa,8.333333333333334)-
        14.784*t10*t11*t12/pow(rhoa,5.666666666666667))+0.444444444444444*
        t1*t15*t2/pow(rhoa,1.666666666666667)-0.180444444444444*t1*
        t2*pow(t20,3.0)*t22*t36+0.28*t1*t2*t20*t22*t28*t32-0.4*t1*
        t16*t2*t21*t32+0.373333333333333*t1*t16*t2*t27*t28-0.133333333333333*
        t1*t2*t20*t21*t26);
    res[6] = 0.5*(-0.180444444444444*t1*t2*t22*t25*t27*t36+
        0.186666666666667*t1*t2*t20*t22*t28*t33-0.266666666666667*
        t1*t16*t2*t21*t33+0.093333333333333*t1*t2*t22*t25*t28*t32-
        0.1*t1*t2*t21*t22*(6.336*t10*t11*grada*t31+19.7037037037037*
        t6*t7*t24*t30+0.0375*t3*t23*t29)+0.248888888888889*t1*t16*
        t2*t20*t25*t28-0.044444444444444*t1*t2*t21*t25*t26);
    res[7] = 0.5*(-0.180444444444444*t1*t2*t20*t22*t34*t36+
        0.093333333333333*t1*t2*t20*t22*t28*t35-0.133333333333333*
        t1*t16*t2*t21*t35+0.124444444444444*t1*t16*t2*t28*t34+0.186666666666667*
        t1*t2*t22*t25*t28*t33-0.1*t1*(-1.728*t10*t11*t19-9.333333333333332*
        t6*t7*t12*t18-0.020833333333333*t3*t8*t17)*t2*t21*t22);
    res[8] = -0.05*t1*t2*t21*t22*(3.5*t6*t7*grada*t9+0.010416666666667*
        t3*t24*t5)-0.090222222222222*t1*t2*t22*pow(t25,3.0)*t36+0.14*
        t1*t2*t22*t25*t28*t35;

}

static void
pw86x_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    pw86x_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw86x_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
pw86x_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/pow(3.141592653589793,4.0);
    t4 = pow(grada,6.0);
    t5 = 1/pow(rhoa,8.0);
    t6 = 1/t1;
    t7 = 1/pow(3.141592653589793,2.666666666666667);
    t8 = pow(grada,4.0);
    t9 = 1/pow(rhoa,5.333333333333333);
    t10 = 1/pow(6.0,0.666666666666667);
    t11 = 1/pow(3.141592653589793,1.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,2.666666666666667);
    t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
        t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
    t15 = pow(t14,0.066666666666667);
    t16 = pow(rhoa,0.333333333333333);
    t17 = 1/pow(rhoa,9.0);
    t18 = 1/pow(rhoa,6.333333333333333);
    t19 = 1/pow(rhoa,3.666666666666667);
    t20 = -0.864*t10*t11*t12*t19-0.777777777777778*t6*t7*
        t8*t18-6.944444444444445E-4*t3*t4*t17;
    t21 = 1/pow(t14,0.933333333333333);
    t22 = pow(rhoa,1.333333333333333);
    t23 = pow(grada,5.0);
    t24 = pow(grada,3.0);
    t25 = 0.648*t10*t11*grada*t13+0.583333333333333*t6*t7*
        t24*t9+5.208333333333333E-4*t3*t23*t5;
    t26 = 1/pow(rhoa,0.666666666666667);
    t27 = pow(t20,2.0);
    t28 = 1/pow(t14,1.933333333333333);
    t29 = 1/pow(rhoa,10.0);
    t30 = 1/pow(rhoa,7.333333333333333);
    t31 = 1/pow(rhoa,4.666666666666667);
    t32 = 3.168*t10*t11*t12*t31+4.925925925925925*t6*t7*t8*
        t30+0.00625*t3*t4*t29;
    t33 = -1.728*t10*t11*grada*t19-3.111111111111111*t6*t7*
        t24*t18-0.004166666666667*t3*t23*t17;
    t34 = pow(t25,2.0);
    t35 = 0.648*t10*t11*t13+1.75*t6*t7*t12*t9+0.002604166666667*
        t3*t8*t5;
    t36 = 1/pow(rhoa,1.666666666666667);
    t37 = pow(t20,3.0);
    t38 = 1/pow(t14,2.933333333333333);
    t39 = 1/pow(rhoa,11.0);
    t40 = 1/pow(rhoa,8.333333333333334);
    t41 = 1/pow(rhoa,5.666666666666667);
    t42 = -14.784*t10*t11*t12*t41-36.12345679012345*t6*t7*
        t8*t40-0.0625*t3*t4*t39;
    t43 = 6.336*t10*t11*grada*t31+19.7037037037037*t6*t7*
        t24*t30+0.0375*t3*t23*t29;
    t44 = -1.728*t10*t11*t19-9.333333333333332*t6*t7*t12*
        t18-0.020833333333333*t3*t8*t17;
    t45 = pow(t25,3.0);
    t46 = 3.5*t6*t7*grada*t9+0.010416666666667*t3*t24*t5;
    t47 = 1/
        pow(t14,3.933333333333333);

   /* code */
    res[0] = 0.5*(-0.1*t1*t2*t20*t21*t22-2.0*t1*t15*t16*t2);
    res[1] = -0.05*t1*t2*t21*t22*t25;
    res[2] = 0.5*(-0.1*t1*t2*t21*t22*t32+0.093333333333333*
        t1*t2*t22*t27*t28-0.666666666666667*t1*t15*t2*t26-0.266666666666667*
        t1*t16*t2*t20*t21);
    res[3] = 0.5*(-0.1*t1*t2*t21*t22*t33+0.093333333333333*
        t1*t2*t20*t22*t25*t28-0.133333333333333*t1*t16*t2*t21*t25);
    res[4] = 0.046666666666667*t1*t2*t22*t28*t34-0.05*t1*
        t2*t21*t22*t35;
    res[5] = 0.5*(-0.1*t1*t2*t21*t22*t42-0.180444444444444*
        t1*t2*t22*t37*t38+0.444444444444444*t1*t15*t2*t36+0.28*t1*
        t2*t20*t22*t28*t32-0.4*t1*t16*t2*t21*t32+0.373333333333333*
        t1*t16*t2*t27*t28-0.133333333333333*t1*t2*t20*t21*t26);
    res[6] = 0.5*(-0.1*t1*t2*t21*t22*t43-0.180444444444444*
        t1*t2*t22*t25*t27*t38+0.186666666666667*t1*t2*t20*t22*t28*
        t33-0.266666666666667*t1*t16*t2*t21*t33+0.093333333333333*
        t1*t2*t22*t25*t28*t32+0.248888888888889*t1*t16*t2*t20*t25*
        t28-0.044444444444444*t1*t2*t21*t25*t26);
    res[7] = 0.5*(-0.1*t1*t2*t21*t22*t44-0.180444444444444*
        t1*t2*t20*t22*t34*t38+0.093333333333333*t1*t2*t20*t22*t28*
        t35-0.133333333333333*t1*t16*t2*t21*t35+0.124444444444444*
        t1*t16*t2*t28*t34+0.186666666666667*t1*t2*t22*t25*t28*t33);
    res[8] = -0.05*t1*t2*t21*t22*t46-0.090222222222222*t1*
        t2*t22*t38*t45+0.14*t1*t2*t22*t25*t28*t35;
    res[9] = 0.5*(-0.1*t1*t2*t21*t22*(0.6875*t3*t4/pow(rhoa,
        12.0)+301.0288065843621*t6*t7*t8/pow(rhoa,9.333333333333334)+
        83.776*t10*t11*t12/pow(rhoa,6.666666666666667))+0.529303703703704*
        t1*t2*pow(t20,4.0)*t22*t47+0.373333333333333*t1*t2*t20*t22*
        t28*t42-0.533333333333333*t1*t16*t2*t21*t42-0.96237037037037*
        t1*t16*t2*t37*t38-1.082666666666667*t1*t2*t22*t27*t32*t38+
        0.118518518518519*t1*t2*t20*t21*t36+0.28*t1*t2*t22*t28*pow(t32,
        2.0)+1.493333333333333*t1*t16*t2*t20*t28*t32-0.266666666666667*
        t1*t2*t21*t26*t32+0.248888888888889*t1*t2*t26*t27*t28-0.740740740740741*
        t1*t13*t15*t2);
    res[10] = 0.5*(0.529303703703704*t1*t2*t22*t25*t37*t47+
        0.28*t1*t2*t20*t22*t28*t43-0.4*t1*t16*t2*t21*t43+0.093333333333333*
        t1*t2*t22*t25*t28*t42-0.1*t1*t2*t21*t22*(-29.568*t10*t11*grada*
        t41-144.4938271604938*t6*t7*t24*t40-0.375*t3*t23*t39)-0.541333333333333*
        t1*t2*t22*t27*t33*t38-0.541333333333333*t1*t2*t20*t22*t25*
        t32*t38-0.721777777777778*t1*t16*t2*t25*t27*t38+0.02962962962963*
        t1*t2*t21*t25*t36+0.28*t1*t2*t22*t28*t32*t33+0.746666666666667*
        t1*t16*t2*t20*t28*t33-0.133333333333333*t1*t2*t21*t26*t33+
        0.373333333333333*t1*t16*t2*t25*t28*t32+0.124444444444444*
        t1*t2*t20*t25*t26*t28);
    res[11] = 0.5*(0.529303703703704*t1*t2*t22*t27*t34*t47+
        0.186666666666667*t1*t2*t20*t22*t28*t44-0.266666666666667*
        t1*t16*t2*t21*t44+0.186666666666667*t1*t2*t22*t25*t28*t43-
        0.180444444444444*t1*t2*t22*t27*t35*t38-0.180444444444444*
        t1*t2*t22*t32*t34*t38-0.481185185185185*t1*t16*t2*t20*t34*
        t38-0.721777777777778*t1*t2*t20*t22*t25*t33*t38+0.093333333333333*
        t1*t2*t22*t28*t32*t35+0.248888888888889*t1*t16*t2*t20*t28*
        t35-0.044444444444444*t1*t2*t21*t26*t35+0.041481481481481*
        t1*t2*t26*t28*t34+0.186666666666667*t1*t2*t22*t28*pow(t33,
        2.0)+0.497777777777778*t1*t16*t2*t25*t28*t33-0.1*t1*t2*t21*
        t22*(6.336*t10*t11*t31+59.1111111111111*t6*t7*t12*t30+0.1875*
        t3*t8*t29));
    res[12] = 0.5*(0.529303703703704*t1*t2*t20*t22*t45*t47+
        0.093333333333333*t1*t2*t20*t22*t28*t46-0.133333333333333*
        t1*t16*t2*t21*t46-0.240592592592593*t1*t16*t2*t38*t45+0.28*
        t1*t2*t22*t25*t28*t44-0.541333333333333*t1*t2*t20*t22*t25*
        t35*t38-0.541333333333333*t1*t2*t22*t33*t34*t38+0.28*t1*t2*
        t22*t28*t33*t35+0.373333333333333*t1*t16*t2*t25*t28*t35-0.1*
        t1*(-18.66666666666666*t6*t7*grada*t18-0.083333333333333*t3*
        t24*t17)*t2*t21*t22);
    res[13] = -0.05*t1*t2*t21*t22*(3.5*t6*t7*t9+0.03125*t3*
        t12*t5)+0.264651851851852*t1*t2*t22*pow(t25,4.0)*t47+0.186666666666667*
        t1*t2*t22*t25*t28*t46-0.541333333333333*t1*t2*t22*t34*t35*
        t38+0.14*t1*t2*t22*t28*pow(t35,2.0);

}

static void
pw86x_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    pw86x_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw86x_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
